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Abstract 

A two-dimensional crystal on the surface of a sphere experiences elastic stress due to the incompatibility of the crystal axes 
and the curvature. A common mechanism to relax elastic stress is the Asaro-Tiller-Grinfeld (ATG) instability. With a combined 
numerical and analytical approach we demonstrate, that also curvature induced stress in surface crystals can be relaxed by 
the long wave length ATG instability. The numerical results are obtained using a surface phase-field crystal (PFC) model, 
from which we determine the characteristic wave numbers of the ATG instability for various surface coverages corresponding 
to different curvature induced compressions. The results are compared with an analytic expression for the characteristic wave 
number, obtained from a continuum approach which accounts for hexagonal crystals and intrinsic PFC symmetries. We find 
our numerical results in accordance with the analytical predictions. 

PACS numbers: 81.10.Aj, 83.10.Rs, 68.08.De 


I. INTRODUCTION 

In material sciences nano structures are of crucial im¬ 
portance, as they often define the macroscopic properties 
of the material. The kinetic effects occurring upon the 
formation of these structures are widely studied, well un¬ 
derstood and often used to control the formation process. 
However, under elastic stress also an interface at rest 
can develop an instability and lead to the formation of 
nano structures. This stress-driven instability, known as 
Asaro-Tiller-Grinfeld (ATG) instability was first studied 
by Asaro and Tiller pQ and later independently by Grin- 
feld [2] and Srolovitz [3]. These authors studied the lin¬ 
ear instability of a planar interface of a stressed solid and 
found that the surface is unstable for perturbations with 
wave numbers less than a critical value. The instabil¬ 
ity is manifested by mass transport. The elastic stress in 
the solid is a destabilizing factor, while the interfacial en¬ 
ergy is a stabilizing one and their interplay leads to inter¬ 
face modulations relaxing the elastic stress. It has been 
extensively studied theoretically and numerically for a 
wide range of different stress driven rearrangement in¬ 
stabilities, see e.g. mm. More recently, the connection 
between the original continuum formulation and a crys¬ 
tal of discrete constituents was successfully established 

[IMS]. 

Less understood is the role of elastic stress, which 
arises from curvature effects for crystals on curved sur¬ 
faces. Such a situation can be found e.g. in the cases of 
coatings, the assembly of biomembranes, the formation 
of molecular monolayers or in the packing of filament 
bundles fTbfH8] . The natural lattice packing of these two 
dimensional crystals is incompatible with the curvature 
of the surface, since the symmetry axes of the crystals are 
bent by the curvature, leading to stressed crystals. The 
influence of this curvature induced elastic stress and its 
relaxation is under investigation in this letter. In m an 
elastic instability of a growing colloidal crystal is consid¬ 
ered experimentally on a spherical droplet and the behav¬ 
ior is analyzed using a continuum theory. Here, we will 



FIG. 1. The temporal evolution (left to right) of a ribbon of 
stressed atoms on the surface of a sphere. The curvature in¬ 
duced stress enhances initially small crystal interface pertur¬ 
bations, which grow exponentially in time. For visualization 
purposes the maxima in the particle density field are extracted 
and considered as atoms. 


instead consider a two-dimensional crystal on a spheri¬ 
cal surface at rest and account for discrete constituents 
of the crystal by using the phase-field crystal (PFC) ap¬ 
proach introduced in nano, see also the review [22] for 
the wide applicability of the modeling approach in hard 
and soft matter systems. The PFC model was also suc¬ 
cessfully applied to crystals on curved surfaces [23h26j . 
but mainly focusing on defects describing grain boundary 
scars m and pleats 28], or properties of Pickering emul¬ 
sions and Bijels [211. A comprehensive investigation of 
curvature induced stress relaxation using the PFC model 
is missing and will be provided in this letter. 

We will briefly introduce the PFC model and the nu¬ 
merical approach and numerically investigate the relax¬ 
ation of curvature induced stress. The results are com¬ 
pared with an analytical continuum model taking into 
account the hexagonal structure of the crystal and the in¬ 
trinsic symmetry of the PFC approach. The comparison 
allows us to identify the relaxation as an ATG instability. 


II. PHASE-FIELD CRYSTAL MODEL 

The phase-field crystal (PFC) model is described by an 
energy functional in terms of the reduced particle density 
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m = j -|VrVf + ]|A r ^| 2 + /(VO dr, (1) 

where /(VO = |(1 — r)jp 2 + \^p A . T denotes the curved 
surface and Vr and Ar are the corresponding surface 
gradient and surface Laplacian. Within our numerical 
consideration T will be a sphere. The only parameters 
are r, corresponding to an undercooling and the average 
density of the system ip. Depending on the parameter 
r the energy functional is minimized by periodic and/or 
constant solutions, modelling a crystal and its melt, re¬ 
spectively m • The temporal evolution of the system is 
given by the H~ x gradient flow 


dt'i/j = Ar 


5F[ip] 

Sip 


( 2 ) 


making ip a conserved quantity. In [30] a derivation of 
this equation from a surface dynamic density functional 
theory (DDFT) is sketched, following the detailed deriva¬ 
tion of the PFC equation in flat space in [2TJI32]. 

Within a one mode approximation in flat space m 
the periodic solution is given in the (x, y )-plain by 


Vv = A 


cos (qx) cos 




( 3 ) 


with equilibrium wave number q and amplitude A defined 
as 


<1 = // A = | (ip + ] /-15r - 36V> 2 ^ • (4) 

Correspondingly, the equilibrium wave length is 
ao = 27 x/q. We will use this one mode approximation 
to initially set up a curvature induced stressed configu¬ 
ration. 


III. NUMERICAL SIMULATION 

Now, we numerically investigate the relaxation of such 
a curvature induced stressed configuration. As in the 
classical ATG instability, only wavelengths of interface 
perturbations above some critical value can be expected 
to grow exponentially. We therefore need large, stressed 
crystals and start with a one mode approximation of 
a ribbon wrapped completely around the equator of a 
sphere with radius R = 100ao/27r (see fig. [l] for an 
overview). Its initial width is adjusted such that it closely 
corresponds to the surface coverage as determined by the 
PFC parameter ip. This allows to keep crystal growth 
at bay in order to avoid competing dynamic instabili¬ 
ties ( e.g . Mullins-Sekerka). As ip steers the width of 
the crystal ribbon and each new particle layer is ex¬ 
posed to increasing stress due to the curvature of the 


sphere, the width can be used to realize setups with dif¬ 
ferently stressed crystals. We choose the PFC param¬ 
eter r = —0.25 to ensure short simulation times and 
'ip = -0.32, 'ip = -0.31, V 5 = -0.30 and ip = -0.29 for 
increasingly stressed crystals. The PFC equations (J2| 
are solved by using a basis decomposition into spherical 
harmonics combined with a semi implicit Euler time dis¬ 
cretization mm- In order to shorten simulation times, 
small amplitude noise was added in each time step. 

Upon temporal evolution, we determine the upper and 
lower interface of the crystal, denoted by b(x) where x 
is a longitude coordinate. The interface mean values are 
constant after an initial relaxation of the initial condi¬ 
tion and before crystalline defects are incorporated at a 
late stage of the interface modulation. The interfacial 
Fourier components b(k) are calculated and the ampli¬ 
tude for each wave number k is monitored over time. 
Subsequently, in order to obtain the growth rate cr k for 
each wave number &, we fit the obtained data to an ex¬ 
ponential function oc exp (a k t) using the time interval of 
constant mean interface. 

We introduce the compression c m = (ao — a) /ao with 
the actual and equilibrium particle distance at the in¬ 
terface a and ao, respectively. In fig. [2] the numerically 
obtained growth rate and interface spectrum b(k) are ex- 
emplarily shown for ip = —0.29. Defining the mean value 
of the interface as the crystal interface, this corresponds 
to a compression at the crystal interface of c m = 18.4%. 
We extract a maximum growth rate of a k = 3.0 x 10 -2 for 
the wave number k ma x = 0.18. The slightly noisy struc¬ 
ture of the growth rates and spectra in fig. [2] originates 
from the noise imposed in the numerical simulations. Un¬ 
fortunately, we also observe considerably large values for 
the growth rate for wave numbers, where the according 
spectra of the interface lacks contributions. These non¬ 
physical contributions are present only due to the rigor¬ 
ous application of our oc e akt fitting procedure, although 
spectral components exhibit no exponential behavior at 
all in this spectral regime. Consequently, these parts of 
the curves should not be taken into account, as suggested 
by the shaded areas in fig. [2] In fig. [ 3 ] we demonstrate 
for the same parameter setting, that the absolute values 
of the maximum growth rate (Jk,max and wave number of 
maximum growth k max do not depend on the number of 
particle layers necessary to realize a certain coverage of 
the sphere (and thus a certain curvature induced stress 
at the crystal interface), as a doubling of the sphere ra¬ 
dius (and thus doubling the number of layers) does not 
considerably alter the previously found values. 

Further data of k max for different compression rates 
Cm appear as solid dots in fig. [5] For the additional data 
c m = 7.5% (ip = -0.32), Cm = 10.2% (ip = -0.31) 
and Cm = 15.3% ('ip = —0.30), we observe increasing 
maximum values of cr k = 1.4 x 10 -4 , a k = 1.8 x 10 -3 
and <Jk = 1-6 x 10 -2 and increasing wave numbers of 
maximum growth rate k ma x = 0.04, k ma x = 0.09 and 
kmax 0.14. 

These increasing growth rates a k and wave numbers 
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FIG. 2. Exemplary interface spectra (a) and growth rates (b) averaged over 20 single noise realizations for a compression 
Cm = 18.4% (ijj = —0.29). The shaded area indicate spectral domains, where the values of the growth rate are non-physical. 



wave number k 


FIG. 3. The equality of growth rates for surface radii R = 
100 ao and R = 200 ao demonstrates their independence on 
the number of particle layers = —0.29). 


instability for an isotropic continuous medium was ana¬ 
lyzed in flat space. 

We start with the general expression for the elastic 
energy (see, e.g., [35] ) 

F = F 0 T ^ Cijkl^ijtkl 5 

where we use the elastic constants Cijki and the strain 
tensors e^- = \{dui/dxj + duj/dxi) with the displace¬ 
ment fields Ui. The indices obey i,j,k,l G {x,y} for two 
spatial coordinates x\ = x, X 2 = y , see fig. [4] We use 
Voigt’s notation xx 1, yy 2 and xy = yx 3. 
Exploiting the intrinsic symmetries for the elastic con¬ 
stants and strain tensors and additionally accounting for 
hexagonal crystal and intrinsic PFC symmetry results in 

Cn = C22 = 3C33 = 3Ci2, Cis = C23 = 0. 


kmax for increasing stress are in accordance with the orig¬ 
inal continuum ATG theory m, as well as with previ¬ 
ous observations from numerical simulations within the 
amplitude equation approach m in flat space. In [14] 
a ’perfect relaxation’ condition was formulated. That 
condition assumes, that a crystal of discrete constituents 
reaches a completely stress free state, when it includes 
a certain number of defects. Equally distributing these 
defects along the crystal interface defines a wave num¬ 
ber. This wave number of maximum stress relaxation is 
plotted as function of the compression c m in fig. [5] and 
nicely agrees with our numerical results, even if the ori¬ 
gin of stress is different. Additionally, we calculate the 
most unstable wave number k ma x within a continuum 
elasticity model. 


IV. CONTINUUM ELASTICITY 

Using the continuum elasticity theory with hexagonal 
symmetry and intrinsic PFC symmetries, we derive ex¬ 
pressions for the wave number of a maximum growth rate 
kmax • We thereby closely follow [34], where the ATG 


Thus, the free energy reads 

F = + (3e^ + 3e 2 yy + Ae 2 xy + 2e xx e yy ) . 

The stress tensor ciij = dF/deij = Cijki^ki is related to 
the strain tensor via 


&ij — 4C33 6ij T ( 2 €ij T dkk) 5 

and vice versa 




^ ~ 4 C 33 + 4C33 U * ^ 
and the equilibrium equation to solve reads 


1 


Sin 


&kk \ 

2 ) 5 


dx k 


which is satisfied for 


(X xx — 


&X 

dy 


2 ’ 


= 0 


dxdy ’ 


d 2 x 


ayy ~ dx 2 


( 5 ) 


with an arbitrary Airy stress function X m x( x iV)- At 
the crystal interface we formulate 


^nn — Pl-> &nt — 0? 


(6) 
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FIG. 4. The geometry of the continuum elasticity model. Un¬ 
der external stress modulations ah ^ of the unperturbed 
interface h ^ grow. 


where pi is the liquid pressure and 

^nn — 'H'i&ij'Tlji &tt — ti&ijtji &nt — W'i&ijtj 

are the normal and tangential components of the stress 
tensor. Describing the crystal interface as y = h(x,t), 
the normal and tangent vectors of the interface are given 

by 


(t x ,t y ) = {l,dh/dx)/yjl + ( dh/dxji , 

( n x ,n y ) = (-dh/dx, l)/\/l + ( dh/dx ) 2 . 

We now solve eq. ([ 6 ]) in a perturbative manner. We 
assume ft to be a small parameter, which describes 
the strength of the interface modulation. The interface 
h(x,t) now reads: 

h = /i (0) + ah w + a 2 /i (2) + • • • . (7) 


Similarly, also the remaining variables are expanded in a 
power series in a 


(o) 

- • • —a . .' 




+ at-- + a €■■ + • • 


2 ( 2 ) 


r (°) 




=^'+a<rg ) +a 2 di ) + - 

Xij =xf + a x\f + a 2 Xif + ■ 


Plugging in the expansions and reordering all terms by 
powers of <a, eq. © reads up to first order 

0 =W + 4s + “ (4» “ 2/i (1), 4°j, ) ) (8) 

o + «(-d* h(iy + °xy + < 4 ° )/l(1) ) > ( 9 ) 

where ' denotes the derivative w.r.t x. Evaluating these 
equations to zeroth order in a gives 

a yJ = ~Pi’ °i° )= 0 (10) 


and (Txx is the applied stress. To proceed with the first 
order in <a, we make the ansatz 

= (yf + By) exp (ikx + ky + c ot) (11) 


and determine the constants A , B by evaluating the in¬ 
terface conditions © and © to first order in a. Assum¬ 
ing further, that the first order perturbation of the flat 
(h^°\x) = const) interface obeys 

h& = h\\ exp ( ikx + uot) , ( 12 ) 


we end up with 

= — 2kaohn exp ( ikx + c ot ), 
a2y =ik(T 0 hn exp ( ikx + ut ), 

=0, (13) 

where we introduced <t (i = <JxJ — ay°J. 

The temporal evolution of the surface perturbation 
h (b is induced by the solidification of liquid at the crystal 
interface. The solidification is driven by the difference of 
the chemical potential between the liquid and solid phase 
A/i Illiquid, l^solid- 


dh dhdh W 

— =-h a - 

dt dt dt 


^=/A„, 


(14) 


with some proportionality constant /. We encounter the 
same thermodynamic situation at the crystal interface 
that is described in detail in |34], appendix A. However, 
we consider hexagonal crystals, making the mathematical 
expressions slightly more extensive. 

We start with the two phases at equilibrium. When 
transforming a small mass element at the interface of 
volume SV from liquid into solid, the change of the Gibbs 
free energy is 


AG = AF + A (piSV) = SV A/i, (15) 


where A F is the Helmholtz free energy change. The 
Helmholtz free energy A F = A Fi + A F m is composed 
of the free energy change of the transformed mass ele¬ 
ment A F m and an interface contribution A F^. The con¬ 
tribution A Fi = 7 kSV accounts for the change of the 
interface free energy caused by the interface tension 7 
and the interface curvature n. For the interface eq. 0 
the curvature k is given to first order a 


hT 

(1 + h ,2 ) 3/2 


—k 2 ah^ 1 \ 


giving 


A Fi = -7 k 2 ah^5V. (16) 

The free energy change of the transformed mass element 
A F m is given by the work necessary to increase its in¬ 
ternal strain to the value of the surrounding solid. The 
work for an infinitesimal change in strain is 

d / = Vijdeij. (17) 
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Together with the change of the volume element d (SV) 
due to increasing strain, the infinitesimal elastic contri¬ 
bution to AG reads 


dG ei =(iijdeij5V Pid(SV ) 

— \{&nnd^nn T O’ttd^-tt) ^nn (^nn T G^tt)] SV 

— ip~tt ^nn) de tt SV. (18) 

We used the mechanical equilibrium eq. © by substitut¬ 
ing —pi = G nn and exploiting a nt = 0, i.e. the stress 
tensor is diagonal in the coordinate system defined by 
the vectors normal and tangential to the interface. We 
express de u and (a tt — a nn ) in terms of the small param¬ 
eter a and discard all terms of higher order than one: 

_ d a tt d ai°J da$ dd 

£ “ ~~ 4C 33 + 8 C 33 8 C 33 + “ 8 C 33 “ 8 C 33 


and 


0 tt - cr ran =<J 0 (l - 2kah ( ' 1 \x,t)'j . (19) 

Using the relations (13) for , d Gyy and integrating 
eq. (18) over stress values up to the interface values given 
in eqs. (10), (13) results in 


A Gel — 


-8C33 


7 — kah 


8C33 


After plugging this together with eq. (16) into eq. (15) 
and dividing by SV, we arrive at 


Afi = 


8C33 


(l-4kahW) + Y ^-(l-2kahW) 


+ 


^0 

8C33 


kah^^j —'yk 2 ah( 1 \ 


( 20 ) 


Comparing to eq. (14) we deduce the exponential growth 
rate uo of the interface perturbation 


w =/ 


(17-U (2i) 


which is maximal for the wave number 
k -- 

r^max — ^ 


8 7 C 33 


( 22 ) 


Now, the result for the wave number of maximum growth 
is adapted to the crystal on the curved surface of a sphere. 


V. DISCUSSION 


In the derivation of eq. (22), all defining quantities were 


evaluated at the crystal interface. In particular, only the 
elastic properties and the chemical potential difference of 


latitude d [rad] 

0.00 0.32 0.45 0.55 0.64 



misfit c m 


FIG. 5. The wave number kmax of the maximum growth 
rate is plotted as a function of the compression c m . For the 
crystals on the sphere, the compression is curvature induced 
and related to the latitude $ of the crystal interface c m = 
1 — cos('d). The results from numerical simulations (solid blue 
line with dots) are bounded by the two limiting cases of zero 
stress (red dashed line) and zero strain (green dashed line) 
in y-direction. The blue dash dotted line corresponds to the 
’perfect relaxation’ condition from m and agrees very well 
with our results. 


the interface determine the value of k max . Thus, we iden¬ 
tify the crystal interface from the previous calculations 
in flat geometry with the interface of the crystal on the 
curved surface of the sphere. Because the elastic constant 
C33 can be obtained directly from the PFC parameters r 
and t/S EQ, the remaining task is to determine the inter¬ 
face tension 7 and the externally applied stress gq with 
regard to the curvature of the spherical surface. We start 
with the strain in x-direction (parallel to the unperturbed 
interface h^) for the crystal on the sphere. The strain 
is e£x = 1 — cos('d) for a crystal interface at latitude $, 
provided the crystal is stress free at the equator (# = 0 ). 
This emphasizes the difference to flat geometry: the cur¬ 
vature induced compression c m = 1 — cos('d) increases 
for larger latitudes, whereas the compression is constant 
in the flat geometry. For the y-direction (perpendicu¬ 
lar to h (°)), we can either assume zero strain (e^ = 0 ) 
or zero stress (< 7 ® = 0). The two cases correspond to 
pi = —c m Cs 3 or to zero liquid pressure pi = 0 , respec¬ 
tively. The interface tension 7 is obtained as the ratio 
of the energy d F needed to prolong the interface by a 
certain length d L and dL, i.e. 7 = dF/dL. We obtain 


a 0 =7 = -CssCm zero stress, (23) 

2 

do =-7 = 2C^c rn zero strain. (24) 

o 

Accordingly, we finally get 

kmax =c m zero stress, (25) 

kmax = \ c -m zero strain. (26) 
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These two lines are also plotted in fig. [5] The actual liq¬ 
uid pressure lies between the two limiting cases of no 
liquid pressure pi = 0 and the case where the liquid 
pressure is so strong that it ensures zero displacement 
in y-direction pi = —c m Cs 3 - Our numerical results are in 
between these two limiting lines and thus allows identify¬ 
ing the observed elastic instability as a curvature induced 
ATG instability. 

VI. SUMMARY 

In summary, we numerically simulated the relaxation 
of curvature induced elastic stress for crystals on a spher¬ 
ical surface within a PFC model. In agreement with 
an analytical continuum model accounting for hexagonal 
crystal and inherent PFC symmetries, the relaxation is 
mediated by the ATG instability. Accordingly, we found 


that the elastic stress at the crystal interface is defin¬ 
ing the growth rates and the characteristic wavelength 
of maximum growth of the interface modulations. This 
situation is different to the elastic instability discussed in 
m, where the growth of a crystal under the influence of 
curvature induced stress leads to anisotropic growth. 

Even if the used numerical approach is restricted to the 
geometry of a sphere using one of the other numerical ap¬ 
proaches discussed in m any curved surface can be con¬ 
sidered or even surface modulations, possibly induced by 
the stressed crystal [36 or resulting from external forces. 
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